# M. Orsi - Essex group - October 2006
# This program opens the file "area_volume.dat", which is one of BRAHMS's output files, and extract the volume compressibility  modulus

# USAGE: cat area_volume.dat | python computeVolumeCompressibilityModulus.py

import sys, string
from math import sqrt
lines = sys.stdin.readlines()

#constants
BOLTZMANN = 1.3806505e-23 # [J / K]

# variables
T = 303.15 # [K]
nLipids = 128  # number of lipid molecules in the system
nWaters = 4232 # number of water molecules in the system
waterVolume = 0.03 # nm^3

print "REMEMBER TO ADJUST THE NUMBER OF LIPID, WATER MOLECULES AND TEMPERATURE!!!"

lineCounter = 0
meanLipidVolumeSum = 0
squareLipidVolumeSum = 0
meanBoxVolumeSum = 0
squareBoxVolumeSum = 0

brahmsAreaVolFile = open( 'area_volume.dat', 'r' )

for line in lines:
    lineCounter = lineCounter + 1
    words = string.split( line )
    boxVolume = string.atof(words[1]);
    meanBoxVolumeSum = meanBoxVolumeSum + boxVolume;
    squareBoxVolumeSum = squareBoxVolumeSum + boxVolume**2
    
brahmsAreaVolFile.close

# computing statistical quantities
meanBoxVolume = meanBoxVolumeSum / lineCounter;
print "meanBoxVolume = %4.3f nm^3" % (meanBoxVolume)
meanSquaredBoxVolumeFluctuation = squareBoxVolumeSum / lineCounter - meanBoxVolume**2
# print "meanSquaredVolumeFluctuation = %6.3f A^4" % (meanSquaredVolumeFluctuation * 10000)

meanLipidVolume = ( meanBoxVolume - waterVolume * nWaters ) / nLipids;
print "meanLipidVolume = %4.3f nm^3" % (meanLipidVolume)

# computing modulus
KV_box = BOLTZMANN * T * meanBoxVolume / meanSquaredBoxVolumeFluctuation  # [ J / nm^3 ]

# conversion considering that J/nm^3 =  pow(10,19) kbar, being given that 1 kbar = 10^8 N/m^2 
KV_box = KV_box * pow(10,19)

print "KV = %6.3f kbar" % ( KV_box )
